CLIMATE 405: Machine Learning for Earth and Environmental Sciences; FALL 2024¶
Mohammed Ombadi (ombadi@umich.edu)¶
Lecture 10 (Wednesday, 10/02/2024)¶
Topics covered in this lecture:¶
- An example of applying Multilayer Perceptron (MLP) for rainfall-runoff modeling.
- Note: Review PowerPoint slides for lecture 10 (Gradient Descent, Backpropagation, Stochastic Gradient Descent, Architechture of a Neural Network, Activation Functions
Import libraries¶
# Import libraries
import os
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.graph_objects as go
import plotly.express as px
import warnings
from scipy.integrate import odeint
from IPython.display import Image
from statistics import mode
from scipy.stats import pearsonr, spearmanr
from sklearn.feature_selection import mutual_info_regression
import statsmodels.api as sm
import pingouin as pg
# Suppress warnings
warnings.filterwarnings('ignore')
# Set number of decimals for np print options
np.set_printoptions(precision=3)
# Set the current working directory
os.chdir(sys.path[0])
Hsu et al. (1995) Study: MLP for Rainfall-Runoff Modelling ¶
This study is the earliest example of applying Multilayer Perceptron (MLP) for rainfall-runoff modelling. The study can be found here
display(Image(filename = 'Hsu et al_WRR 1995.png', width= 600, height= 300))
Objectives:
- The objective of this study is to use a multilayer perceptron (MLP) to predict (simulate) runoff using rainfall data as an input.
- Compare the performance of the MLP model to an autoregressive moving average with exogenous variables (ARMAX)
- Compre the performance to a physics-based conceptual model; The Sacramento Soil Moisture Accounting Model (SAC-SMA)
A diagram desribing the SAC-SMA model:
display(Image(filename = 'SAC-SMA.png', width= 600, height= 300))
Case Study:
- The Leaf River Basin, near Collins in Mississipi (size = 1949 km2)
- Calibration data: 1983 water year (10/01/1982 - 09/30/1983)
- Validation data: water years 1984 - 1988
The General Architecture of the MLP model:
display(Image(filename = 'Fig 1_ Hsu et al_WRR 1995.png', width= 600, height= 300))
Splitting the data into training and testing:
display(Image(filename = 'Fig 2_ Hsu et al_WRR 1995.png', width= 600, height= 300))
Model Selection:
The inputs consist of:
- Lagged Rainfall(t-$\tau$), $\tau = 1, 2, 3, 4, 5$
- Lagged Runoff(t-$\tau$), $\tau = 1, 2, 3, 4, 5$
In the table below:
- ANN(5, 4, 3, 1) means a multilayer perceptron with 5 lagged inputs for rainfall, 4 lagged inputs for runoff, hidden layer with 3 nodes, and 1 output.
display(Image(filename = 'Table 1_ Hsu et al_WRR 1995.png', width= 800, height= 600))
Validation:
display(Image(filename = 'Fig 3b_ Hsu et al_WRR 1995.png', width= 800, height= 600))
Example: Building an MLP for Rainfall-Runoff Modelling ¶
Let's build an MLP model to replicate the Hsu et al (1995) study. Here, we're going to use data:
- From a watershed in Maine (Fish River Watershed)
- For the time period 01/01/2016 to 12/31/2021
Let's upload and plot the data:
fish_river = pd.read_csv('fish_river.csv')
print(fish_river.shape)
fish_river
(2192, 4)
| datetime | year | rainfall | runoff | |
|---|---|---|---|---|
| 0 | 2016-01-01 | 2016 | 0.000675 | 1720.0 |
| 1 | 2016-01-02 | 2016 | 0.144395 | 1680.0 |
| 2 | 2016-01-03 | 2016 | 2.918217 | 1640.0 |
| 3 | 2016-01-04 | 2016 | 0.000000 | 1580.0 |
| 4 | 2016-01-05 | 2016 | 0.000000 | 1520.0 |
| ... | ... | ... | ... | ... |
| 2187 | 2021-12-27 | 2021 | 0.000000 | 989.0 |
| 2188 | 2021-12-28 | 2021 | 3.711825 | 972.0 |
| 2189 | 2021-12-29 | 2021 | 0.018728 | 947.0 |
| 2190 | 2021-12-30 | 2021 | 0.000000 | 928.0 |
| 2191 | 2021-12-31 | 2021 | 0.030013 | 902.0 |
2192 rows × 4 columns
Processing the data:¶
Let's process the input data to have 4 variables (rainfall lagged by 1 and 2 times steps, and runoff lagged by 1 and 2 time steps):
data = fish_river.copy()
data['rainfall_tau_1'] = data['rainfall'].shift(1)
data['rainfall_tau_2'] = data['rainfall'].shift(2)
data['runoff_tau_1'] = data['runoff'].shift(1)
data['runoff_tau_2'] = data['runoff'].shift(2)
data = data.dropna()
data
| datetime | year | rainfall | runoff | rainfall_tau_1 | rainfall_tau_2 | runoff_tau_1 | runoff_tau_2 | |
|---|---|---|---|---|---|---|---|---|
| 2 | 2016-01-03 | 2016 | 2.918217 | 1640.0 | 0.144395 | 0.000675 | 1680.0 | 1720.0 |
| 3 | 2016-01-04 | 2016 | 0.000000 | 1580.0 | 2.918217 | 0.144395 | 1640.0 | 1680.0 |
| 4 | 2016-01-05 | 2016 | 0.000000 | 1520.0 | 0.000000 | 2.918217 | 1580.0 | 1640.0 |
| 5 | 2016-01-06 | 2016 | 0.000000 | 1450.0 | 0.000000 | 0.000000 | 1520.0 | 1580.0 |
| 6 | 2016-01-07 | 2016 | 0.000000 | 1390.0 | 0.000000 | 0.000000 | 1450.0 | 1520.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2187 | 2021-12-27 | 2021 | 0.000000 | 989.0 | 0.000000 | 0.000000 | 1020.0 | 1070.0 |
| 2188 | 2021-12-28 | 2021 | 3.711825 | 972.0 | 0.000000 | 0.000000 | 989.0 | 1020.0 |
| 2189 | 2021-12-29 | 2021 | 0.018728 | 947.0 | 3.711825 | 0.000000 | 972.0 | 989.0 |
| 2190 | 2021-12-30 | 2021 | 0.000000 | 928.0 | 0.018728 | 3.711825 | 947.0 | 972.0 |
| 2191 | 2021-12-31 | 2021 | 0.030013 | 902.0 | 0.000000 | 0.018728 | 928.0 | 947.0 |
2190 rows × 8 columns
# Inputs (4 variables)
X = data[['rainfall_tau_1', 'rainfall_tau_2', 'runoff_tau_1', 'runoff_tau_2']]
X
| rainfall_tau_1 | rainfall_tau_2 | runoff_tau_1 | runoff_tau_2 | |
|---|---|---|---|---|
| 2 | 0.144395 | 0.000675 | 1680.0 | 1720.0 |
| 3 | 2.918217 | 0.144395 | 1640.0 | 1680.0 |
| 4 | 0.000000 | 2.918217 | 1580.0 | 1640.0 |
| 5 | 0.000000 | 0.000000 | 1520.0 | 1580.0 |
| 6 | 0.000000 | 0.000000 | 1450.0 | 1520.0 |
| ... | ... | ... | ... | ... |
| 2187 | 0.000000 | 0.000000 | 1020.0 | 1070.0 |
| 2188 | 0.000000 | 0.000000 | 989.0 | 1020.0 |
| 2189 | 3.711825 | 0.000000 | 972.0 | 989.0 |
| 2190 | 0.018728 | 3.711825 | 947.0 | 972.0 |
| 2191 | 0.000000 | 0.018728 | 928.0 | 947.0 |
2190 rows × 4 columns
# Ouputs (1 variable)
y = data[['runoff']]
y
| runoff | |
|---|---|
| 2 | 1640.0 |
| 3 | 1580.0 |
| 4 | 1520.0 |
| 5 | 1450.0 |
| 6 | 1390.0 |
| ... | ... |
| 2187 | 989.0 |
| 2188 | 972.0 |
| 2189 | 947.0 |
| 2190 | 928.0 |
| 2191 | 902.0 |
2190 rows × 1 columns
Building the model:¶
We're going to build a multilayer perceptron neural network for prediction, using Scikit-learn MLPRregressor. The documentation of the function is here
from sklearn.neural_network import MLPRegressor
from sklearn.model_selection import train_test_split
# Create an MLP (one hidden layer with 5 nodes)
mlp = MLPRegressor(hidden_layer_sizes=(5,),
activation= 'relu',
solver= 'adam',
learning_rate= 'constant',
learning_rate_init= 0.001,
alpha= 0,
max_iter=1000,
random_state=42)
# Split the data into training and test sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
# Fit the model to the training data
mlp.fit(X_train, y_train)
MLPRegressor(alpha=0, hidden_layer_sizes=(5,), max_iter=1000, random_state=42)In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
MLPRegressor(alpha=0, hidden_layer_sizes=(5,), max_iter=1000, random_state=42)
Evaluating the model:¶
pred = mlp.predict(X_test)
obs = y_test
obs = obs.values.reshape(-1,)
plt.plot(np.arange(len(obs)), obs, 'or', np.arange(len(pred)), pred, 'k')
[<matplotlib.lines.Line2D at 0x31ffc9ac0>, <matplotlib.lines.Line2D at 0x31ffca750>]
plt.scatter(obs, pred)
<matplotlib.collections.PathCollection at 0x31f14f7a0>
print('RMSE = ', np.round(np.sqrt(np.sum( (obs-pred)**2))))
print('Pearson correlation = ', np.round(np.corrcoef(obs, pred)[0,1], 3))
print('Relative Bias = ', np.round((np.mean(pred-obs)/np.mean(obs))*100), 3, ' %')
RMSE = 5956.0 Pearson correlation = 0.993 Relative Bias = -0.0 3 %
Extract Training Information:¶
Loss Curve
Loss_curve = mlp.loss_curve_
print(len(Loss_curve))
plt.plot(Loss_curve)
1000
[<matplotlib.lines.Line2D at 0x3200e7800>]
Number of samples for each iteration
Activation Functions ¶
def sigmoid(x):
y = 1 / (1 + np.exp(-x))
return y
x = np.arange(-10,11)
y = sigmoid(x)
plt.plot(x, y)
plt.grid()
plt.savefig('sigmoid.jpeg', dpi= 600)
def ReLU(x):
x = pd.Series(x)
y = x.apply(lambda x: max(0, x))
return y
x = np.arange(-10,11)
y = ReLU(x)
plt.plot(x, y)
plt.grid()
plt.savefig('ReLU.jpeg', dpi= 600)
def softmax(x):
y = np.exp(x) / np.sum(np.exp(x))
return y
x = np.random.random_integers(1, 14, 7)
y = softmax(x)
plt.subplot(1,2,1)
sns.barplot(x)
plt.grid()
plt.subplot(1,2,2)
sns.barplot(y)
plt.grid()
plt.savefig('softmax.jpeg', dpi= 600)
plt.show()